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DEVEL(»’I€NT OF A DOSIICTER FOR DISTRIBIH^D BODY ORGANS 


This report includes the progress of research in the study of develofasMit 
of a dosimeter for distributed body organs. 

The basis for the great interest in the development of a space dosimeter 
is the dose which a particle will deposit in human tissue. In the attached 
paper, the calculational methods for estimation of dose from external proton 
exposure of arbitrary convex bodies is briefly reviewed and all of the necessary 
information for the estimation of dose in soft tissue is presented. The effects 
of nuclear reaction which become important for determining dose equivalent are 
included in these calculations. This work on "Proton-Tissue Dose Calculations" 
is proposed for publication as a NASA-TO. 

The above results are currently being applied towards the development of 
space radiation dosimeti/ of distributed body organs. 



Proton-Tissue Dose Calculatims 

John W. Wilson ^ 

NASA- Langley Research Center 
Hampton, Va. 

and 

Govind S. Khandelwal 
Old Dominion University 
Norfolk, Va. 

ABSTRACT 

Calculational methods for estimation of dose from external proton 
esqposure of aribtrary convex bodies is briefly reviewed and all of the* 
necessary information for the estimation of dose in soft tissue is 
presented. Special en^hasis is on retaining the effects of nuclear 
reaction especially in relation to the dose equivalent. 

INTRODUCTION 

IVhen an object is exposed to external radiation, the dose field. 
within the object is a complicated function of the character of the 
external radiation, the shape of the object (including orientation), 
and the object’s material con^osition. Calculation of dose i^ithin an 
object involves solution of the appropriate Boltzmann transport equatitm 
where the external radiation source imposes boundary conditions on the 
solution. Althou^ general purpose computer programs exist for making 
such estimates (ref. 1), they are seldom used in practice when the object 
is bounded by a complicated surface as, for example, is the human body. 
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Instead, calculatlcms are usually Bade for single gemaetric shapes 
froB whidi inferences are then nade for sore general geometries and 
the resultant errors are uncertain. 

In the case of external proton radiation such as that encountered 
near hi^«energy accelerators, in space, and in high-altitu^ air- 
craft, it was found that the problem of dose estimation could be greatly 
simplified (ref. 2) and still include the effects of nuclear reaction, 
which in^ses the major hurdles in any accurate calculation, with a hi^ 
degree of accuracy. Furtherrore, it was shown that the method, tAen in 
error, was always conservative. Required for such calciilations is a 
Icnowledge of the transition of protons in semi-infinite slab geometxjr 
which is the simplest geometry for existing transport computer prograas-. 
Indeed, almost everything that is known about the dose in, humans to 
external proton radiation is inferred from calculations in slab 
geometry (ref. 3). 

In the present note, a general method for estimation of dose in 
arbitrary convex geometry in teJms of dose conversion factors in 'slab, 
geometry is briefly discussed. These dose conversion factors for protons 
in tissue are then represented using buildtq) factors. A parametric form 
for the buildup factors is presented , *.hThe values for the parameters 
are derived from Monte Carlo calculations of various authors. All of 
the necessary information to estimate dose and dose equivalent for proton 
irradiation of convex objects of arbitrary shape is contained heroin. 
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THEORY 


In passing throu|^ tissue, energetic protons interact mostly 
throu^ ionization of atomic constituents by the transfer of small 
amounts of momentum to orbital electrons. Althou^ the nuclear 
reactions are far less numerous, their effects are magnified because 
of the large momentum transferred to the nuclear particles and the 
struck nucleus Itself. Unlike the secondary electrons formed throuih 
atomic ionization by interaction- with the primary protons, the resulting 
radiations of nuclear reaction are mostly heavily ionizing and generally 
have large biological effectiveness. Many of the secondary particles 
of nuclear reactions are sufficiently energetic to promote similar 
nuclear reactions and thus cause a buildup of secondary radiations. 

The description of such processes requires solution of the transport 
equation. The approximate solution for the transition of protons in 
30 cm thick slabs of soft tissue for fixed incident energeis are pre- 
sented in references 4 through 11* The results of such calculations 
are dose conversion factors for relating the primary monoenergetic 
proton fluence to dose or dose equivalent as a function of position in 
a tissue slab. 

Whenever the radiation is spatially uniform, the dose at any point 
in a convex object may be calculated according to reference 2. by 

.OCX ) s* ^ f (Ij 

where R^(z,E) is the dose at depth z for normal incident protons of 
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energy E on a tissue slab« ^(n«E) is a differential proton flusnce 
along direction and is the distance from the boundary 

along iT to the point It has been shown that equation (1) always 
overestimates the dose, but is an accurate estimate when the ratio 
of the proton beam divergence due to nuclear reaction to the bodies 
radius of curvature is small. Equation (1) is a practical prescription 
for introducing nuclear reactim effects into calculations of dose in 
geometrically complex objects as the human body. The main requirearat 
is that the dose conversion factors for a tissue slab be adequately 
known for a broad range of energies and depths* 

Available information on conversion factors are for discrete 
energies from 100 MeV to 1 TeV in rather broad energy steps and for 
depths from 0 to 30 cm in semi-infinite slabs of tissue (zefs. 4,5,8, 
and 9) . The nuclear reaction data used for high'Onergy nuclecms is 
usually based on Monte Carlo estimates (refs. 12-14) with low-energy 
neutron reaction data taken from esqjerimental observation. The quality 
factor as defined by the ICRP (ref, 15) is used for protons. The 
quality factor for heavier fragments and the recoiling nuclei is axbi- 
trarily set to 20 which is considered conservative although the average 
quality factor obtained by calculation is comparable to estimates 
obtained through observations made in nuclear emulsion (ref. 16). 

To fully utilize equation (1), the fluence-to-dose conversion factors 
for normal incident protons on a tissue slab must be known for all 
energies and depths. A parametri ration of the conversion factors was 
introduced by Wilson and Khandelwal (ref. 2) which allowed reliable 


hmcEsuo »Aae slamk tror mtm 
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interpolation and extrapolation fron known values. In tiie follouing^ 
a refinement aiul extension of that work will be discussed* 

« ‘ I » I. vi; »i* 

Fluence-tO'Dose Conversion Factors 


Ihe conversion factor I^(s»E) is composed of two terns represen- 
ting dose due to tiie primary bean protons* md tiie dose due to secondary 
parades produced in nuclear reaction. *. IbuSt 





where the pxinary dose equivalent conversion factor is giysn by 


&) = ?(B) / P(e*J 

The LET denoted by SCE) in equation (3) is calculated using Bathe's 
formula above 243.8 keV as given by 


where 


S(S) = 


Hmt/o e* Z 




(4a) 


2 ■. average atomic number 
A ■ average atomic wei^t 
I - adjusted ionization potential 
m B electron nass 
e B electron charge 
V B proton velocity 
c B velocity of light 
N^b Avogadro's number 

At proton energies below 243.8 KeV, the LET is calculated by the 

enpirical esqtression 
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(4b) 


which approximately accounts for the inner shell corrections in toft 


tissue. The proton range in soft tissue is given by 

c 

H(B)= \ d*'/S<0 

with the reduced energy in equation (3) given by 

"a] • 


(S) 


( 6 ) 


where €(>:) is inverse function of R(E). The total nuclear survival 
probability for a proton of energy E is given by 


where the macroscopic cross section o(E) for tissue as calculated by 
Bertini is given by Alsmiller et al. (ref. 18). The proton total optical 
thickness given by 

<rce')JcVs<e9 (8) 

is tabulated in table 1 for purposes of numerical interpolation. &i 
the case of conversion factors for absorbed dose, the R„(z,E) is 
taken as 

R/i,S)-=. ?(£) S(^)/F(e^) 




(») 
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Bulldii|> Pacton 

The representation of the conversion factors is sin^lified (see 
ref. 2) by rewriting equatim (2) as 

idiere P(z«E) is recognized as the dose buildi^ factor. Hie aaia 
advantage for introducing the buildi? factor into equation (10) is that 
unlike R^(z,E), the buildiq) factor is a saoothly varying function, oi: 
energy at all depths in the slab and can be aj^roximated by the sisple 
functim. 

F(i, s) = ( A, + ^ e + 

idiere the parameters are understood to be energy dependent. Ihe 
coefficients are found by fitting equation (11) to the valiws of 
the buildup factors as estimated from the Monte Carlo calculations of 
proton conversion factors. The resulting coefficients are shown in 
table 2. The coefficients for 100, 200, and 300 MeV protons were 
obtained using the Monte Carlo data of Thmer et al. (ref. 4). the 
values at 400 , 730, 1500 and 3000 MeV were obtained from Hie results 
of Alsmiller and Armstrong (ref. 9). The 10 GeV entry was obtained 
from the calculations of Armstrong and Chandler (ref. 9). Values noted 
in table 2 by asterisk on the corresponding energy were obtained by 
interpolating between data points or smoothly extrapolating to imit 
buildup factor at proton energies near the Coulcmib barrier for tissue 
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nuclei (- 12 MeV). The resulting buildup factors are shown in figures 
1 and 2 in coi^aiison to the MoniiJ Carlo results where the error ban 
were deternined by drawing smooth limiting curves so as to bracket the 
^k)nte Carlo values and to follow the general functional ^pendMice. 

These uncertainty limits should, therfore, be interpreted as approximately 
Zo limits, rather than lo ranges usually used in expressing tsicertainty 
limits. 


CONVERSION FACTOR COMPUTER COOi 

To utilize equation (1) in a specific problem requires values 
for the conversion factor R^(z,E) over the range of interest. Formtla# 
for these factors are presented in the previous section. A cc»q>utev 
code has been generated to return values of R^(z,E) for arbitrary depth 
z and energy E. This code is listed in the appendix and is described 
briefly here. There are six main functions to be generated relating to 
LET, range- energy relations, quality factor, and the functions relating 
to nuclear reaction effects given as nuclear survival probability and 
buildup factor. 

The functions relating to ionization by the primary beam are 
generated by the function subroutine RTISS. Tables for R(E), and S(E) 
are generated on the first call to RTISS. Subsequent intermediate 
values are found by numerical interpolation above 10 KeV. A sicplified 
approximation based on equation (4b) is used at lower energies. The 
function c(x) is found by numerical inversion of R(C). 
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1)16 quality factor is approximated by 

OpCS) St (12) 

for S gt^ater than 3S ^V/cm and set to unity for smaller LET. 

The values shovm in table 1 of the c^tical itonsity are generated 
in ^e functim subroutine PN(E) and stored in an array for numerical 
interpolation and the nuclear survival probability is calculated using 
equation (7). . 

The coefficients for calculating the buildup factors axe generated 
by subroutine ANTER as a function of energy by interpolating between 
the valu&£ t^hown in table 2. 

Ihc cctiversion factors are generated by subroutine RESP by supply 
parameters z and £ which represent distance in centimeters of 
tissue and proton energy E in units of MeV. The returned values of 
the conversion factors have units of rad (or rem], per proton per 
centimeter squared. 
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SAMPLE CALCULATlCm 


To illustrate the usage of the buildup factors described here« 
calculations of the dose in slab geometry for normal incident .protons 
with spectra typical of the space environment have beat made. Cal- 
culations were also made neglecting nuclear reaction effects and the 

t. 

percentage contribution to the dose and dose equivalent due to 
nuclear reactions are shown in figures 3 and 4. The spectra indicated 
by GCR in the figures represent galactic cosmic radiation with 
spectrum given by 

'ecu 


^ -f-0/fit^f ( 13 ) 

The spectra denoted by the parameter represent solar cosmic ray 
spectra given ais * 




(14) 


with the rigidity given as 


Pf^) = (15) 

where q is the proton charge and is the proton mass. The value 
Pq “ 100 MV corresponds to an intermediate-energy solar event :ahd 
P^ « *.00 MV corresponds to a high-energy solar event. The curve denoted 
by » 100 MeV represents the energetic inner belt protons with spectrum 


(piE) ( 1 ?) 

It is clear from the figures that dose estimates for galactic cosmic 
rays and high energy solar cosmic rays cannot be accurately calculated 
without proper account of nuclear reactions. Ihis is especially trt» 
for estimat'; of the dose equivalent. 


Although reasonable estimates (± 10%) of low and intermediate 

solar cosmic ray absorbed doses are e:q>ected, the dose equivalent 

estimates must include nuclear reaction effects. Marginally good 

estimates of absorbed dose for inner belt protons can be made by 

neglecting nuclear reactions but dose equivalent estimates require 

inclusion. nuclear reaction effects. 

A 
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APPENDIX 


PROGRAM LISTING FOR CONVERSION FACTOR CALCULATION 


FRBClSDim HIAHK NOT lUiOO 


SUBROUTINE RESPCEN*X.RAD*REMJ 

C THIS SUBROUTINE GENERATES VALUES FOR THE SLA3 CONVERSION FACTORS 

C FOR values of PROTON ENERGY EN <MEV) AND DEPTH IN THE SLAB X CCM) 

real C(8) 

ENER*EN . 

EXaX 

CALL ANTERIENER*C»RJ 

RRSS®R-EX 

ENERPsETlSiRRES) 

IF(ENERP)34*33*34 
33 continue 

RAQsO* 

REMaO* 

RETURN 

3A continue 

CALL APR0B(EX*ENER*PR0B) 
call AT0PPCENERP*ST0PP| 

2 call AF(ST0PP»0ALF) 

22 PES=PR08*ST0PP#0ALF 

C 0 R£Q=(C {1 )-»-X*{C(2)+X*CC3) n*EXP(-X*CC4)| 

COREC® (C(5>+X*{C(6)+X*-C(7) ) )*EXPC-X*CC8) ) 

IF (COREQ.lt. I • ) COREO»I. 

IF(cOREC.LT.l . ) COREC*l. 

REM=PES*CCREO »1.6E~8 
PESsOR08*ST0PP 
RAD=PES»C0REC*1 .6E-8 
RETURN 

end 


SUBROUTINE AnTER (ENER*C.R) 

C THIS subroutine GENERATES THE VALUES OF THE PARAMETERS 

C OF THE ANALYTIC FITS OF THE MONTE CARLO RESULTS 

real C(8).A(12«8)«E(l2l 
logical FALS 

data E/30.t60..100.. ISO. « 200 • .300 • » 400 • . 730. * 1200*. 1500*. 3000 •• 
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I lOOOO*/ 

DATA A/I .O. 1 *2*1 »A« 1 .5*1 *6* 1 *70* 1 .90 • 3* AO • 4 • 32 «4*63 «5»35« 6»20« 
2 0»0*0»0»*C2*»07»*09»*lI»*I3»*15f>«*I67«»lT0*»l 90 » • 2Q0 * 

3C.0, .Ot 0*0. 0.0,0. 0»0.C*-7.0» . 0003^> , • CO 1 AS . .0025 . .0030 , .DOaS* 
AO.C, «C1 3. .030. .0385* .0 AO , .033 « .0228 , .01 50 • • 0 1 3 * .0 1 2 . .0 1 0 • .01 0. 
®*i.0.|.0«l.l.l*l2*l.|5.l.2*l.24«i.4*1.6V«l.S.2*.2#w. 

69 .0. .0 1 « .040, .06. •062. .063* .071 , *o9. .094 « .095. .1 0. • 1 1 * 

70. 0. 0. 0«0.0«0.0 *0.0*0. 0.0. O. .00^1 , ,00080* .0015* »0C 2 *.00205. 
BO.C* .01 ,.026* .031 ,.032, .026* .0228,. 01 5 *.0122, .01 2*. 01. .01/ 

DATA FALS/.T./ 
data IPT/-1/, 

«sPTIS(ENE»l 
IF(FAUS) GO TO 10 
1 CONTINUE 

EL0S«AL06(ENEP) 

CALL lUNl < 12* 12.E.8* A.2.ELOG*C. IPT. lERR) 

RETURN 

10 CONTINUE 

no 11 1*1.12 
E(n*AL06(E( 1 M 

11 CONTINUE 
FALS*»F» 

GO TO 1 
END 


SU8R0UTINE AF (STOPP.QALF) 

C this subroutine computes the quality FACTOR AS A FUNCTION OF 

C LINEAR ENERGY TRANSFER 

IF(ST0PP-35.)1 1.11 .12 

11 0ALF*1. 

RETURN 

12 0ALF=.06*ST0PPA*.8 
RETURN 

END 


SUBROUTINE APROB (EX .E .PROS I 

C THIS subroutine generates values for the NUCLEAR SURVIVAL PROBABILITY 


OF A PROTON OF ENERGY E CMEV> AFTER TRAVELING A DISTANCE EX (CM) IN TISSUE 


SRES»RT 1S(E)-EX 
PROB-0* « 

IF (RRES.LE#0, ) RETURN 
ENEW-ETISCRRES) 

PROS =PN ( E ) /PN { ENEW ) 

RETURN 

FNH 


FUNCTION PN(E) 

PN GIVES PROBABILITY THAT PROTON TRAVELS FULL RANGE WITHOUT 
BEING ABSORBED 

external fox 

LOGICAL TRU 
REAL R(30)*ET(3O) 

data ET/O* , 10. *25* *50« * 100* * 150» *200 • • 25U* ♦ 300 • * 350 • «400* «500«« 

1700. . 900. * I 100.. 1300.. 1500. , 1 700 . .2000 . .2200 .. 2400 . .2600. «2800«* 

23000. . 4000. *5000. .6000. .7000.. 3500. . 10000./ 

DATA TRU/.T./ 

DATA 1PT/-1/ 

IF(TRU) GO to lO 
111 FR*E 

call IUNI <30.30 .ET. 1 .R. 2. ER. BYRD. IPT. lERR) 

PNs EXP (-BYRD) 

RETURN 
10 TRUs.F. 

RtllsO. 

DO 1 1=2.30 
EU=ET(I) 

GsET( I-l ) 

CALL MGAUSS(G.EU.04.ANS.FOX.F. I ) 

R{ I )=R( I-l )+ANS 
1 CONT I Nue 
PRINT 19 

19 FORMAT (///.25X.*PN GRID*//) 

PRINT 119 

119 FORMAT (10X.*E VALUES FOR GRID*//) 

PRINT 226. ET 
226 FORMAT (2X.8E15.6) 




PPINT 227 

227 format t//*lOX«»R values FOR G*?I0*) 
PRINT 226«R 
GO TO III 
END 


SUSROUTINE FOX(X«P) 
ENERsX 

3 CALL ASIGMIENER. SIGMA! 
call ATOPP{ENER*STOPP» 
F*S|GMA/STOPP 
2 RETURN 
END 


SU3R0UTINE ASIGM (ENEP*SIGMAt I 

1 

C THIS SUBROUTINE GENERATES VALUES OF TOTAL NONELASTIC MACROSCOPIC I 

C CROSS SECTION (CM»*2/G) IN TISSUE AS A FUNCTION OF PROTON ENERGY ENERCMEVI | 

real EN ( A3) •CR0SC43) | 

data EN/25«32*29,86»34# 16«39«86*AA.65t50*0j .6Q.19*70*2A*79*A7*89»9 | 

1 I t 100.8t 1 1 7.9» 139. 3t 156.3. 175#3» 186.6. 202. 9»266. 1 »304.7 i 375#2*407. j 

27.471.6.507.1 .574.5.61 1 . 4 . 678.3 . 7l 4 , 5 . 776.4 . 809. 3 . 870.4 . 9 16.8* 1 007 | 

3. . I 129. . 1406* . 1756. .2024. ,2318. .3071 • *2403. .3943. .5000. .8000*. | 

4100C0./ I 

DATA CROS/2.614 .2.360*2. 153.1 .985. 1 .387. 1 .757. 1.621.1 .526.1 .431 .!• | 

1 379. 1 .327, 1.261.1.211.1.187,1.164,1.152.1.141.1 .097. 1 .087.1 .100. t • 

21 36 * 1 . 199. 1 .21 2 . 1 .266* 1 .293. 1 .350. 1 .379. 1 .424. 1 .440. 1 .471 . 1 .478.1 • | 

3504,1 .477. 1 .490 . I .433. 1 .455. 1 .497 , 1 .475, 1.461.1 .463. 1.46.1 .458. | 

4!. 452/ I 

data IPT/-1/ I 

i: e=eneR 

1FIeNER.lt. 25. 32) ENFR=23.32 | 

call IUNI (43, 43*EN. 1 .CR0S.2.ENER. CROSS. IPT. IERR) 

S IGVA*(CROSS/lOO. ) 

ENERsE 

RETURN i 

END 
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FUNCTION RTISCE) 


THIS SU3R0UTINE GENERATES THE RANGE-ENERGY RELATIONS AND LET FOR 
SROTONS IN tissue 

EXTERNAL ATOE 

real ETCS7l*RT<57)tST<57> 

logical false 

data FALSE/.T./ • ■ 

data NP/57/ 

DATA £T/*01 • *02 • • 03 » «04 • *05 • aOS* *07 • *08t *09 «»2*«3*«A«» S* 

1 •6««7f •B* m9 9 1 • «2« *3» 9^9 *5* «8« *9* * 1 0* tEO* «30« ♦ AO« «50 •’t 

260.«70» vaO* «9U*« lOO.tlSO. «200.«300> * A00« • 500* « 600* « 700* * 

3800. *900« *1000 •« 1500. «2000* t25U0« f30C0»«A000*«50C0*»6O0O»* 
47000 .♦ 8500# • I 0000 ./ 

N*l • ‘ . 

IF (FALSE )G0 TO I 0 t 

12 continue. 

RT I S»E »* • 697/ C 25 1 7* ♦ .697 J 

if(e.lt..oi) return 
AaAl_0G(E> • 

DO 1 IE»2*NP ' , . - 

1F(A.LT.ET( IE)) GO TO 2 . 

1 CONTINUE • 

2 i=ie * 

SL0PE*(RT( I )-RT( l-I ) )/(ET( I )-ETC l-l ) ) * 

RAL=RT( I-l )+SLOPE*( A-ET( I-l ) ) « 

RTIS*EXP(RAL) 

RETURN . 

ENTRY STIS . 

N»2 • 

IF (FALSE )G0 TO 10 . * 

13 CONTINUE * 

RTIS*E*».303* (25I7.-6283.*E) 

IF(E*LT..01) RETURN 
AaALOG(E) • 

DO 3 IE=2.NP 

IF( A.LT.ET( IE) ) GO TO 4 . 

3 CONTINUE ♦ 

4 1*16 • 

SLOPE* (ST ( I )-ST ( I-I ) >/(ET( I )-ET( I-l ) ) « 

sal *ST( I-l )+SLOPE»(A-ET(I-l )) • 



0 ^ in 


»TtS*EXPf5ALI 
RETURN • 

entry ETIS 

N«3 f 

lFCrALSE)G3 TO 10 • 

continue t 

RTIS«C2517#».697*E)»*1 .43472 
IF(E.UT..On RETURN 
R«At_OG(E) • 

DO 5 1R>2«NP 

IFCR.LT.RTI IR> ) 60 TO 6 * 

CONTINUE • 

I*IR • 

SL0PE»<ET(I J-ETd-l ) >/(RT(I l-RTC l-l 11 « 

EAL *ET( I-I l + SLOPE*IR-RTCI-n) t 

RTIS-EXPIEALI 
RETURN • 

10 CONTINUE • 

RTCl )«0. 

ST(1 )*0* 

,Vs06 

CO 21 I»2.NP 

call ATOPPCETCll.STlin 

call MGAUS3(ET( I-l 1 .ETII ).M*ANS.AT0E«F*1 I 
21 RTU )«RT( I-l l+ANS 
RIRST*RT<2l 
EIRSTs*ET(2l 
00 U IX»2.NP 
ET( IX)«AL00CET( IXI 1 • 

RTC IX)aALOG(RT< IXI) • 

11 ST ( IX)»ALOG(ST{ IX) 1 • 

false* .F. • 

GO TO (12.13tl4)N • 

ENO • 





SUBROUTINE AT0E(E*F1 

call AT0PP(E*S> 

F*|./S 

return 

ENO 
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S'J«1R0UTINE ATOPP(ENER«STOPP) 

C THIS SU3R0UTINE COMPUTES THE STOPPING POWER FOP PROTON IN TISSUE 

{F(EN£P*CT**243S) GO TO 2 

STOPP> (251 7« -6283 •*EnEP )*ENEP*<^» 303 

RETURN 

2 2ETA*ENER/938*2ll 

3ETaS*( CZETA*CZETA-*-2* I )/C (ZETA-M « )**2U 
V/SEal .022201£6»8ETAS/(1 t-BETAS) 

FSET aALOG ( WBE ) -BET AS 

5TOPPa>30726lAB»(-2*2378342-»>«529726*F6ET)/BETAS 

RETURN 

END 


SUBROUT I NE MGAUSS ( A • 0 • N . SUM • FUNC ♦ FOFX ♦ NUMBER > 


DIMENSION U(5 J tR(5> tSUMC I t »FOFXC 1 ) 
00 I UL«1. NUMBER 
I SUM(LU)*0#0 

1F(A.EQ*B) RETURN 
U( 1 la«425562830S09|B4 
U ( 2 }»• 283302302935376 
U(3)a« 1 6029521 S8504B8 
U ( 4 ) a • 0 674683 1 6655508 
U < 5 > ■ • 0 1 304673574 I 4 1 4 
R(1 )s# 1477621 12357376 j 
R ( 2 ) ■ • I 3463335965499 
R(3)a* 109543181257991 
R (4 )»#074 725674575290 
R ( 5 ) « •033335672 1 54344 
FINE«N 

DELT A»F I NE/ ( 0-A ) 

00 3 K»1 «N 
XI ix-i 

fine»a+xi/oelta 
00 2 I 1*1 t5 
UU*U( I I i/delta+fine 
call FUNC (UU«FOFX) 


00 2 JOYSOY*! *NUySER 

? SUM(^OYBOY)«R( I I ) »FOrX( JQY90Y )^SUM( J0Y80Y ) 
00 3 JJ«| «S 

UU*{ 1»0-U< JJ» I/OELTA+FINE 
call FUNC <UU*F0FXI 
DO 3 NN« It number 
3 SUM ( NN ) «R f J J ) *FOFX ( NN ) 4SUM C NN ) 

00 7 fJK«ltNUMBER 
7 SUM ( UK )aSUM(lJK)/OELTA 
RETURN 
END 
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SW3*50UTINC I UN I US£S MRST OR SECOND ORDER 
LftORANGfAN INTCR' wLAT I ON TO ESfliMATt T^-£ VALUES 
0- A SET Or functions AT A »0|MT XO* IWNI 
USES ONE INOESFNOENT VARIABLE TAfiLE. ANO A OEPRNOENT 
7AR1A9LE TASlE FOR EACH FUNCT I ON TO 3E KVALUATEO# 

the routine accepts the independent variables spaced 
AT equal or unequal intervals* each dependent* 
variable table must contain function values corres- 
ponding TO EACH XCll IN the JKDEPENOENT VARIABLE 
TAIlF* the estimated values are returned in the VO 
ARRAY WITH THE N-TH VALUE OF THE ARRAY HOLDING THE 

value of the n-th Function value evaluated at xo* 

• A ’ ' “ ‘ * 

call IUNI «NMAX*N**X*faTAB#Y* IOROER*3iCQ*YO*IPT* lERR) 


F I UN 1 003 
♦IUNIC04* 
^ I UN i cost 
♦IUNI006C 
♦lUNlOCTt 

♦ itwiooat 

*IUN1009< 

♦lUNlOlOf 

FlONtOl tc 

*1UNI012C 

*lUNtOl3C 

♦lUNlOlAC' 

♦lUNlOiSC 

*IUNI0I6C 

♦IUNI017C 

FiUNio'isd 

F-IUNI0I90 

» I UN I 0200 

*tUNt02IO 


I f 

I I 

c* 

^ I c* 
c* 
c* 
c* 
f c* 

■ c* 
^ c» 
: c* 
c* 
: c* 
I c» 
: C* 

C» 

c* 

' c* 

c* 

c* 


CF 

:* 

:* 



* 

* 

|» 

» 



NMAX 


N 


X 


NTA8 

Y 


{ORDER 


XO 


VO 


THE MAXIMUM NUMBER OF POINTS IN THE INDEPENDENT 

VARIABLE ARRAY* 

THE ACTUAL NUMBER OF POINTS IN THE INDEPENDENT 

ARRAY»WHERE N *LE* NMAX* 

A one -dimensional array* dimensioned CNMAXI in THE 

calling program* which contains the independent 

VARIABLES* THESE VALUES MUST BE STRICTLY MONOTONIC. 

THE number Or DEPENDENT VARIABLE TABLES 

A two-dimensional array DIMENSIONED CNMAX.NTABI IN 

THE calling program* EACH COLUMN OF THE ARRAY 

CONTAINS A dependent VARIABLE TABLE 

interpolation parameter supplied by the user* 

*C zero order. INTERPOLATION9 THE FIRST FUNCTION 
VALl'E IN EACH DEPENDENT VARIABLE TABLE IS 
ASS.GNEO TO the CORRESPONDING MEMBER CF THE YO 
array* the FUNCTIONAL VALUE IS ESTIMATED TO 
REMAIN CONSTANT AND EQUAL TO THE NEAREST KNOWN 
FUNCTION VALUE* 


FIUNIC22C 
*1 UN 10230 
*1UNI02AC 
♦I UN I 0250 
» I UN 1 0260 
♦IUNI 0270 
*1 UN I 0280 
♦I UN I 0290 
*1 UN 10300 
♦IUNI0310 
•* I UN 1 0320 
♦1 UN I 0330 
♦I UN I 0340 
♦I UN I 0350 
♦I UN I 0360 
♦IUNI 0370 

♦ I UN I 0380 
♦I UN I 0390 
F| UN 10400 
♦IUNI0410 
♦I UN I 0420 

* I UN 10430 
•I UN I 0440 
♦1 UN 1 0450 
♦I UN 1 0460 


the INPUT POINT AT WHICH I NTERPO-AT ION WILL BE 

per-ormeo* 

A CNE-0 1 MENS I ONAL ARRAY DIMENSIONED CNTABi IN TH£ 
calling PEOGRAM, upon RETURN THE ARRAY CCNTaInS 
ESTIMATED value OF EACH FUNCTION AT XO* 


♦I UN I 0470 
F IUNI 0480 
FI UN I ,.*4 70 
F I'JN I OS “ 0 
FlUNlOSlO 
the FI UN I 0520 
♦IUNT053C 
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IPT 0.N ”-^£ FIPST call IPT MUST 3£ INlTfALlZFO TO -1 SO ^lUMlOSSO 
TH^T VONOTONICI TV WILL BE Ch£C<EO* UPON LEAVING TH£ ^^lUNlC560 
routine IPT EQUALS THE VALUE OF THE INDEX OF THE X *IUNI0570^ 
VALUE PRECEDING XO UNLESS EXTRAPOLATION WAS ♦ I UN 1 0380 

PERFORVEO. in that case the value of IPT IS *IU\I0590 1 

RETURNED AS9 *IUNID500 ' 

=U denotes XC *LT. X(l) IF THE X ARRAY IS IN *IUNI0610= 

increasing order and X(1) •GT* XO IF THE X ARRAY *IUNI0620 ^ 
IS IN DECREASING ORDER. *10X10630 | 

®N DENOTES XO .GT. X(N) IF THE X ARRAY' IS IN *IUN!l05AO 1 

. . increasing order AND XO ,»L,T. X(N) IF .THE X ARRAY *IUN I '3350. I 


IS IN DECREASING ORDER. 

ON SUBSEQUENT CALLS. IPT IS USED AS A POINTER TO 
BEGIN THE SEARCH FOR XO. 

lERR ERROR PARAMETER GENERATED BY THE ROUTITC 

=0 normal return ‘ * 

=J THE J-TH element OF THE X ARRAY IS OUT OF ORDER 
=“1 ZERO ORDER INTERPOLATION PERFORMED BECAUSfe ' 

I ORDER =0. 

= -2 ZERO ORDER I NTERPOLAT I ON PERFORMED BECAUSE ONLY 
ONE POINT WAS. IN X ARRAY. 

=-3 NO interpolation WAS PERFORMED BECAUSE 

INSUFFICIENT POINTS WERE SUPPLIED FOR SECOND 
ORDER INTERPOLATION. 

=-A EXTRAPOLATION WAS PERFORMED 

UPON return THE PARAMETER I ERR SHOULD BE TESTED IN. 

the calling program. 


REQUIRED ROUTINES 


* SOURCE CMP3 ROUTINE MTLUP MODIFIED *IUN 

* BY computer SCIENCES CORPORAT I ON* I 'j\ 

*• *IUN 

* language fortran *IUN 

*■ *IUN 

* * I UN 

* DATE released AUGUST 1.1973 * I UN 

*• #IUN 

* LATEST revision JUNE 9, 1975 *IUN 

* * . *IUN 


NONE 


* lUNl 

*IUNl 

♦lUNl 

*IUN 

* I UN 1 

* I UN 1 
♦tUNi 
♦lUNl 
♦I UNI 
*IUN1 

* I UN 1 
♦lUNl 

* I UN 1 

* I UN I 
*ruN] 

* I UN 1 

* I UN i 

* I UN 1 

* I UN 1 
*IUN1 
♦I UNI 
» I UN 1 

* I UN 1 


language 


date released 


latest revision 


OI '-'ENS ION X( I ) .YINMAX. I ) .YO ( 1 ) 

NMl=N-l 

irpoiiO 

J-l 

TEST “OR ZERO ORDER INTERPOLATION 

DELX=X(2)-X( 1 ) 

I" (I ORDER .EO, 0 ) GO TO lO 
IF Cs.lT. 2) GO rO 20 
C-C T-' 50 








f 


1 j 
20 


30 


40 


50 


C 

c 

c 

.C 

c 


c 

c 

c 


60 


C 

c 

c 


65 


70 

80 

C 

c 

c 


CO Tn 30 

OD 40 NT*1 .\T43 
YO (NT )^Y< 1 .NT ) 

CONTINUE 
5£Tu;?N 

IF CIPT .GT. -!) GO TO 65 

CHECK FOR TABLE OF NODE POINTS BEING STRICTLY MONOTONIC 

the SIGN OF DELX SIGNIFIES WHETHER TABLE IS IN 

INCREASING OR DECREASING ORDER* • 

• . . ■ • . • 

IF COELX .EQ. 0) GO TO 190 - - 

IF (N .EO. 2) GO TO 65 

CHECK for SIGN CONSISTENCY IN THE DIFFERENCES OF 
subsequent PAIRS . .... . . 

DO 60 JsB.NMl . 

IF (DELX * (X( J+1 )-X( J) n 190.190.60 
CONTINUE 

IPT IS INITIALIZED TO BE WITHIN THE INTERVAL 

IF (IPT ,LT. 1) IPT=1 
IF (IPT .GT. N.MI ) IPT=*NM1 
IN= SIGN (I.O.JDELX *( Xy-X(IPT).)I 
o= X( IPT) - XO 

1“ (P* (XdPT +1)- XO)) 90.160.80 
IPT =IPT +IN 

TEST TO SEE IF IT IS NECCESARY TO EXTRAPOLATE 

IF (T T.GT.O .and. IPT .LT. N) GO TO 70 
IERR*-4 

ipt-sipt- in 

TEST -CR order OF I nTERPOLATION 


IF (lOROER .GT, 1) GO TO 120 


FIRST ORDER INTERPOLATION 


lUNl 110^ 

lUNlU t€i 
I UN 1 1 1 2C| 
lUNl I |3ll 
lUNt 1 14C| 
I UN I 1 1 5C| 
lUNll 16CI 
lUNl 1 17C| 
I UNI 11 8C 
lUNl 1 I9d 
IUNII 2 OO 1 
IUNU2lO| 
lUNll^ed 
IUNll230j 
IUNII24Ci 
IUNI1250| 
I UN 112601 
1UNI1270^ 
I UN 11 280^ 
lUNI 1290 
I UN I 1300 
IUNI1310 
I UN I 1320 
IUNI1330 
I UN II 340 
I UN 1 1 350 
I UN I 1360 
IUNI1370 
I UN I 1380 
lUNl 1390 
IUNI1400 
lUNl 1410 
lUNI 1420 
IUNI1430 
lUNl 1440 
lUNI 1450 
IUNI1460 
I UN 11470 
I UN I I 480 
I UN I 1490 
I UN I 1500 
IUNI1510 
I UN 1 1520 
I UN I I 530 


i 


I 


1 

C , ' 



IUNll54f 

1 

t 

IPT1=IPT+1 


lUNi I55f 



XTMPl*XO-X( IPT > 


|UNII56( 

1 


XTMP2=X< IPTl )->; « IP7 ) 


I UNI 157 

j 


XT.ypisxrMPi/xrvaa 


1 UN r 1 58: 

1 

1 


00 100 NTsl.NTAi 


lUNi 159C 

1 

w 


YTMP = Y( :=>”! .\T>-Y( IPT.NT ) 

• « 

lUNI I60f 

u 


YC (NT ) C I = T»N7 )+YTYP»XT.VPl 


iUNll61( 

r 

G- 

too • 

CONTI N'c-E .... . . ... 

• 

tUNIl62f 

1 

1 


IF (lERR .E3. -4| IP7=IPT-MN 


lUNI 163^ 

t 

I 


RETURN 


tUNIl64: 

t 

1 

c 

• 


lUNI I65( 

1' 

c 

SECOND ORDER INTERPOLAT I O-N 


lUNl 166( 

1 

c 



IUNI167C 

1 

120 

IF (N .EO* 2) GO TO 200 


lUNl 168c 

1 

C 



lUNI 169C 

1 

C 

CHOOSING A THIRD POINTSO AS TO MINIMIZE THE 

distance 

lUNl 170C 

1 

c 

between THE THREE POINTS USED TO INTERPOLATE 


IUNI171C 

1 

c 


* 

lUNl 172C 

I 


IF CIPT .EOi NMI ) GO TO 140 


lUNl I73C 

I 


IF (IPT .EO. 1) GO TO 130 


IUNT174C 

1 ' 

— 

A1sABS<X 0-X( lPT-1 ) ) 


lUNi 175C 

1 


AasA3S(X{ IPT+2)-X0) 


IUNI176C 



IF C A 1 -A2 } 1 40 G 1 30 , 1 30 

. • ' 

IUNI177C 


130 

L*IPT 


IUNI17SC 

|- 


GO TO 150 


■ lUNI 179C 

% 

140 

L=IPT -1 

. 

lUNiiaoc 


150 

Vl=X(L)-X0 

■ ; 

. luNiiaic 



V2sX(L+l)-X0 


IUNI182C 



V3=X(L+2)-X0 


lUNl 183c 

t 


DO 160 NT=1,NTAB 


lUNliaAC 

: 


YYls(Y(L.NT) * V2 - Y(L+1*NT) * Vl)/(X(L+1) - X(L)) 


IUNI185C 

S' 


YY2*(Y(L+1 »NT)»V3-Y (L+2tNT ) *V2 } / { X ( Lf2 )-X (L+ 1 ) ) 


IUNI186C 



Y0(NT)s=(YY1»V3-YY2»V1 ) / ( X (L+2 )-X ( L ) ) 


IUNI187C 


160 

continue 


lUNi 188C 

1 


IF (lERR .EQ. -4) IPT=IPT + IN 


lUNI 189C 



RETURN 


lUNI 190C 

w 

ISO 

IF(P *NE. 0) IPTsIPT +1 


1UNI191C 



DO IS5 NTs 1. NT A3 


IUNI192C 



YO(NT)=Y{ IPT«NT 1 


lUNl 1930 


185 

CONTINUE 


lUNI 1940 


RETUf^N lUN 11950 

IUNII960 

IERR is set 70 The subscript of the member of The table lUNri970 


WHICH IS OUT OF ORDER 

is 

^OOR QUALITY 

IUNI19B0 
lUNI 1990 

lERRsJ +1 

IUNI2000 

return 


IUNI2010 

IERRs-3 


IUNI2020 

RE TURN 


lUN 12030 

END 


lUN 12040 






Tissue Optical Thickness for Protons 
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Fifuro 2.- Re* buildiq> factor for several depths in tissue as a function of incident proton energy 
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